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Abstract. We present and analyze a new stable space-time Isogeometric 
Analysis (IgA) method for the numerical solution of parabolic evolution 
equations in fixed and moving spatial computational domains. The dis¬ 
crete bilinear form is elliptic on the IgA space with respect to a discrete 
energy norm. This property together with a corresponding boundedness 
property, consistency and approximation results for the IgA spaces yields 
an a priori discretization error estimate with respect to the discrete norm. 
The theoretical results are confirmed by several numerical experiments 
with low- and high-order IgA spaces. 
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1 Introduction 

Let us consider the parabolic initial-boundary value problem: find u : 
Q —^ M such that 

dtu{x, t) — Au{x, t) = f{x, t) for (x, t) £ Q := f2 x (0, T), 

u{x,t)=0 for (x, t) G := (9l7 X (0, T), (1-1) 

u(x,0) = uo(x) forx G 17, 

as the typical model problem for a linear parabolic evolution equation 
posed in the space-time cylinder Q = 17 x [0,T], where dt denotes the 
partial time derivative, A is the Laplace operator, / is a given source 
function, uq are the given initial data, T is the final time, and 17 C 
{d = 1, 2,3) denotes the spatial computational domain with the boundary 
df2. For the time being, we assume that the domain 17 is fixed, bounded 


2 


U. Langer, S. E. Moore, M. Neumiiller 


and Lipschitz. In many practical, in particular, industrial applications, 
the domain 17, that is also called physical domain, is usually generated 
by some CAD system, i.e., it can be represented by a single patch or 
multiple patches which are images of the parameter domain (0,1)'^ by 
spline or NURBS maps. 

The standard discretization methods in time and space are based on 
time-stepping methods combined with some spatial discretization tech¬ 
nique like the Finite Element Method (FEM) |26I52I31] . The vertical 
method of lines discretizes first in time and then in space [52], whereas 
in the horizontal method of lines, also called Rothe’s method, the dis¬ 
cretization starts with respect to (wrt) the time variable [3T]. The later 
method has some advantages wrt the development of adaptive techniques. 
However, in both approaches, the development of really efficient adaptive 
techniques suffers from the separation of the time and the space dis¬ 
cretizations. Moreover, this separation is even more problematic in paral¬ 
lel computing. The curse of sequentiality of time affects the construction 
of efficient parallel method and their implementation of massively parallel 
computers with several thousands or even hundreds of thousands of cores 
in a very bad way. 

The simplest ideas for space-time solvers are based on time-parallel 
integration techniques for ordinary differential equations that have a long 
history, see m for a comprehensive presentation of this history. The 
most popular parallel time integration method is the parareal method 
that was introduced by Lions, Maday and Turinici in [35]. Time-parallel 
multigrid methods have also a long history. In 1984, Hackbusch proposed 
the so-called parabolic multigrid method that allows a simultaneous exe¬ 
cution on a set of successive time steps m- Lubich and Ostermann [36] 
introduced parallel multigrid wavefrom relaxation methods for parabolic 
problems. A comprehensive presentation of these methods and a survey of 
the references until 1993 can be found in the monograph [56]. Vandewalle 
and Horton investigate the convergence behavior of these time-parallel 
multigrid methods by means of Fourier mode analysis m- Deshpande et 
al. provided a rigorous analysis of time domain parallelism m- Very re¬ 
cently, Gander and Neumiiller have also used the Fourier analysis to con¬ 
struct perfectly scaling parallel space-time multigrid methods for solving 
initial value problems for ordinary differential equations m and initial¬ 
boundary value problems for parabolic PDEs [H]. In these two papers, the 
authors construct stable high-order dG discretizations in time slices. In 
[39] this technique is used to solve the arising linear system of a space-time 
dG discretization, which is also stable in the case of the decomposition 
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of the space-time cylinder into 4d simplices (pentatops) for 3d spatial 
computational domains. This idea opens great opportunities for flexible 
discretizations, adaptivity and the treatment of changing spatial domains 
in time |4()l41l39l25j . We also refer to |55I1I6I49I5()I51I48I37I8I3T] where 
different space-time techniques have been developed. One class of special 
space-time methods, that is very relevant to single-patch space-time IgA, 
which we are going to consider in this paper, is the multiharmonic or 
harmonic balanced FEM that was first used for solving non-linear, time- 
harmonic eddy current problems by electrical engineers, see, e.g., [58]. 
Lately the multiharmonic FEM has been applied to parabolic and eddy 
current time-periodic boundary value problems as well as to the corre¬ 
sponding optimal control problems, see |4l27l28l32j . Babuska and Janik 
already developed h — p versions of the finite element method in space 
with p and h — p approximations in time for parabolic initial bound¬ 
ary value problems in the papers [2] and [3], respectively. In jUj, Schwab 
and Stevenson have recently developed and analyzed space-time adaptive 
wavelet methods for parabolic evolution problems, see also [9]. Similarly, 
Mollet proved uniform stability of an abstract Petrov-Galerkin discretiza¬ 
tions of boundedly invertible operators and applied this result to space- 
time discretizations of linear parabolic problems [38]. Very recently Urban 
and Patera have proved error bounds for reduced basis approximation to 
linear parabolic problems |54] . whereas Steinbach has investigated con¬ 
form space-time finite element approximations to parabolic problems [46j . 
We here also mention the papers by Olshanskii and Reusken who have 
developed Eulerian space-time finite element method for diffusion prob¬ 
lems on evolving surfaces [43142] . Our approach uses special time-upwind 
test functions which are motivated by a space-time streamline diffusion 
method |18l21l22l23j and by a similar approach used in |5|. 


The increasing interest in highly time-parallel space-time methods is 
certainly connected with the fact that parallel computers have rapidly 
developed with respect to number of cores, computation speed, memory, 
availability etc, but also with the complexity of the problems the people 
want to solve. In particular, the optimization of products and processes 
on the basis of computer simulations of the underlying transient processes 
(PDE constrains) foster the development of space-time methods since the 
optimality system is basically nothing but a system of primal and adjoint 
PDEs which are coupled forward and backward in time, see, e.g., [53|. To 
the best of our knowledge, we are not aware of any paper on space-time 
IgA for evolution equations. 
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In this paper, we present a stable discrete space-time variational for¬ 
mulation for the parabolic initial-boundary value problems of the form 
(tnn in the sense that the discrete bilinear form a/i(-, •) : V^h x Vq/i —^ 1^ 
is elliptic on the IgA space Vq/i wrt a special discrete norm \\-\\h- For sim¬ 
plicity, we consider the single-patch case where the space-time cylinder 
Q, that is called physical domain, can be represented by one smooth, 
uniformly regular spline or NURBS map of the parameter domain 
Q = (0,1)'^+^. In IgA, that was introduced by Hughes, Cottrell and 
Bazilevs in 2005 [20], we use the same basis functions for both repre¬ 
senting the approximate solution and defining the geometrical mapping 
Approximation, stability and error estimates for /i-rehned IgA meshes 
of spatial computational domains can be found in [7], see also the mono¬ 
graph |10] for a comprehensive presentation of the IgA and its math¬ 
ematical analysis. Using these approximation results for B-splines and 
NURBS, the Vo/i-ellipticity of the discrete bilinear form a/i(-,-) wrt the 
discrete norm || • ||/i, an appropriate boundedness result and the consis¬ 
tency, we derive asymptotically optimal discretization error estimate in 
the discrete norm || • ||/i. Furthermore, we consider moving spatial compu¬ 
tational domains f2{t) C M'^, t G [0,T], which always lead to a fix space- 
time cylinder Q = {{x,t) G : x G G (0,r)} C that 

can easily be discretized by IgA. We again derive a stable space-time IgA 
scheme and prove asymptotically optimal discretization error estimate in 
a similar discrete norm || • \\h,m- Finally, we present a series of numeri¬ 
cal experiments for hxed and moving spatial computational domains that 
support our theoretical results. Since the discrete bilinear form is Voh- 
elliptic, one may expect that multigrid methods can efficiently solve the 
resulting space-time system of algebraic equation. Indeed, this is the case 
as one of our experiments with the standard AMG code hypre shows in 
Subsection 16.,11 where we solve a linear system with 1.076.890.625 space- 
time unknowns (dofs) within 156 seconds using 16.384 cores. This results 
demonstrates the great potential of our discrete space-time formulation 
for the implementation on massively parallel computers. 

The remainder of the paper is organized as follows: In Section [2l we 
recall the standard space-time variational formulations. Section [3] is de¬ 
voted to the derivation of a stable space-time IgA discretization of our 
parabolic initial-boundary value problem. In Section 0] we derive our a 
priori discretization error estimate. Section [5] deals with the case of mov¬ 
ing spatial computational domains. Our numerical results for fixed and 
moving spatial domains are presented and discussed in Section [6j Finally, 
we draw some conclusions in Section [T] 
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2 Space-Time Variational Formulations 

Let us first introduce the Sobolev spaces = {u £ L 2 {Q) : d^u £ 

-^ 2 (< 3 )) Vet with 0 < |ct| < I, dlu £ L 2 {Q),i = 0, of functions de- 

hned in the space-time cylinder Q, where L 2 {Q) denotes the space of 
square-integrable functions, a = (cti,etd) is a multi-index with non¬ 
negative integers |a| = ai + ... + ad, d^u := d^°‘^u/dx'^ = 

and dlu := d'^u/dt^, see, e.g., [29I30| . Furthermore, 
we need the spaces = {u £ L 2 {Q) : Va;U £ \L 2 {Q)Y‘,u = 

001117} and = {u £ L 2 {Q) : £ [L 2 {Q)]'^, dtu £ L 2 {Q),u = 

0 on S, and tt = 0 on Ut} for introducing the weak space-time formula¬ 
tion of where Ut := f? x {T}, and VxU = {du/dxi,... ,du/dxd)~^ 

denotes the gradient with respect to the spatial variables. 

The standard weak space-time variational formulation of dni reads 
as follows: find u £ such that 

a{u,v) = l{v) Vu £ iLQ’i(Q), (2.1) 

with the bilinear form 

a{u,v) = — J u{x,t)dtv{x,t)dxdt + j Vxu{x,t)-Vxv{x,t)dxdt (2.2) 
Q Q 

and the linear form 

l{v) = J f{x,t)v{x,t)dxdt + J UQ{x)v{x,0)dx. (2.3) 

Q Q 

The variational problem (12.11) is known to have a unique weak solution 
[29130] . In the later books, more general parabolic initial-boundary value 
problems, including other boundary conditions and more general ellip¬ 
tic parts, and non-linear versions are studied. Moreover, beside existence 
and uniqueness results, the reader also hnds useful a priori estimates and 
regularity results, see also m- We here mention that time-stepping meth¬ 
ods rather based on line variational formulations which are formulated in 
function spaces of abstract functions mapping the time interval (0, T) into 
some Sobolev space of functions living on 17, see, e.g., [57159] , [52], and 
[53] for the connection of these different formulations. 

3 Stable Space-Time IgA Discretizations 

In this section, we briefly recall the definitions of B-splines and NURBS 
basis functions and their use for both the geometrical representation of the 
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space-time cylinder Q and the construction of the IgA trial spaces where 
we look for approximate solutions to our variational parabolic evolution 
problem ()2.1j) . For more details on B-splines and NURBS-based IgA, we 
refer to the monograph [TO]. Then we derive our stable space-time IgA 
scheme that uniquely defines an approximate solution Uh in the IgA space. 
This approximate solution can be dehned by solving one linear system of 
algebraic equations the solution of which is nothing but the vector ui^ of 
control points for Uh- 


3.1 B-splines and NURBS 



Fig. 1. Space-time IgA paraphernalia for Q C with d = 1 (left) and 
d = 2 (right). 


Let p > 1 denote the polynomial degree and n the number of basis 
functions defining the B-spline curve. Then a knot vector is nothing but a 
set of non-decreasing sequence of real numbers in the parameter domain 
and can be written as S = {^i,... ,^„+p_|_i}. By convention, we assume 
that ^1 = 0 and ^^-i-p+i = 1- We also note that the knot value may 
be repeated indicating its multiplicity m. We always consider open knot 
vectors, i.e., knot vectors where the first and last knot values appear 
p + 1 times, or, in other word, the multiplicity m of the hrst and last 
knots is p-|- 1. For the one-dimensional parametric domain (0,1), there is 
an underlying mesh KLh consisting of elements K created by the distinct 
knots. We denote the global mesh size by h := max{/ii^ : K E /C/j}, where 
hx '■= diam{K) = length(Ar). For the time being, we assume that the 
ratio of the sizes of neighboring elements is uniformly bounded from above 
and below. In this case, we speak about locally quasi-uniform meshes. 
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The univariate B-spline basis functions Bi^p : (0,1) —)• M are recur¬ 
sively defined by means of the Cox-de Boor recursion formula as follows: 


BiAO 


hpiO 


1 if < ^ < ^i+i, 
0 else, 


Ci+p ~ 




Ci+p+1 ^ 


1(0) 


(3.1) 

(3.2) 


where a division by zero is defined to be zero. We note that a basis 
function of degree p is {p — m) times continuously differentiable across a 
knot with the multiplicity m. If all internal knots have the multiplicity 1, 
then B-splines of degree p are globally continuous. In all numerical 

examples presented in Section El we exactly consider this case. 

Now we define the multivariate B-spline basis functions on the space- 
time parameter domain Q := (0,1)*^^^ C with d = 1,2, or 3 as 

product of the corresponding univariate B-spline basis functions. To do 
this, we first define the knot vectors Hq, = ,... AA+Pa+i} every 

space-time direction ol = 1,... , d -|- 1. Furthermore, we introduce the 
multi-indicies z := (fi,..., id+i) and p := (pi,... ,prf+i), the set Z = {i = 
(zi,..., id+i) : Za = 1, 2,..., Uq,; a = 1, 2,... , d -|- 1}. The product of the 
univariate B-spline basis functions now gives multivariate B-spline basis 
functions 

d-l-l 

Bi,p{0 ■— Bi^^pAA)i ( 3 - 3 ) 

where ^ = (^i,..., .^rf+i) G Q. The univariate and multivariate NURBS 
basis functions are defined in the parametric domain by means of the cor¬ 
responding B-spline basis functions {BiAi^x- givenp = (pi,... ,Pd+i) 
and for all z G X, we define the NURBS basis function Ri^p as follows; 


R,,p;(0,l)^M, (3.4) 

with the weighting function 

1U:(0,1)^M, W{0.= Y.wAA0, (3.5) 

iex 

where the weights Wi are positive real numbers. 

The physical space-time computational domain Q C is now de¬ 
fined from the parametric domain Q = (0, l)'^'*'^ by the geometrical map¬ 
ping 

$ ; Q ^ Q = (k(Q) c M''+\ $(0 = RiAO^i, 

iex 


(3.6) 
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where Ri^p, i £ I, are the multivariate NURBS basis functions and {Pijjgj C 

]^d+i 

are the control points, see also Figure [Das well as Figures [6] and [7] 
in Subsection 16.21 

By means of the geometrical mapping ()3.6p . we define the physi¬ 
cal mesh fCh, associated to the computational domain Q, whose ver¬ 
tices and elements are the images of the vertices and elements of the 
corresponding underlying mesh JCh in the parametric domain Q, i.e., 
JCh := {K = ^{K) : K £ fCh}- We denote the global mesh size of 
the mesh in the physical domain by /i := max{/i;^ : K £ with 

hK '■= \\^^\\Lao{K)hj^ and hj^ = diam(iF). Further, we assume that the 
physical mesh is quasi uniform, i.e., there exists a positive constant Cu, 
independent of h, such that 

hx < h < CufiK for all K £ JCh- (3-7) 

Finally, we dehne the NURBS space 

Vh = span{(^fe,i = Ri^p o $ (3.8) 

on Q by a push-forward of the NURBS basis functions, where we assumed 
that the geometrical mapping <h is invertible a.e. in Q, with smooth inverse 
on each element K of the physical mesh K £ JCh, see [7] and m for more 
details. 

Furthermore, we introduce the subspace 

Voh ■= Vh n hI’^{Q) = {vh £ Vh ■■ VhlEuEo = 0} = span{iph,i}iei (3.9) 

of all functions from Vh fulfilling homogeneous boundary and initial con¬ 
ditions. For simplicity, we below assume the same polynomial degree for 
all coordinate directions, i.e. pa = p for all a = 1, 2,... , d -|- 1. 


3.2 Discrete Space-Time Variational Problem 


In order to derive our discrete scheme for defining an approximate solution 
Uh £ Voh, we multiply our parabolic PDF (II.1|) by a time-upwind test 
function of the form Vh + OhdtVh with an arbitrary Vh £ Uofe U Hqq{Q) 
and a positive constant 9 which will be defined later, and integrate over 
the whole space-time cylinder Q, giving the identity 


+ 0hg,n,8,v, - Au,iv, + OW.v,)) « 

Q 


j fivh + OhdtVh) dxdt 


(3.10) 


Q 
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that is formally valid for all Vh £ Voh C Hq’q{Q). Integration by parts 
with respect to x in the last term of the bilinear form on the left-hand 
side of (I3.10p gives 


Q 

Ux ■ ^xUh {vh + OhdtVh) ds 
OQ 

We mention that dtVh is differentiable wrt x due to the special tensor 
product structure of Voh- Using the facts that Vh and dtVh are always 
zero on E, and the x-components Ux = (ni,..., n^)~^ of the normal n = 
{ni,... ,nd,nd+i)~^ = {nx,nt)~^ is zero on Eq and Et, we arrive at our 
discrete scheme: find Uh G Voh such that 

ah{uh,Vh) = lh{vh) Vu/i G Voh, (3.11) 



^xVh + OhVxUh ■ VxdtVh) dxdt 


j f„ui.At +0k j 

Q Q 


where 


ah{uh,Vh) = j {dtUhVh + OhdtUhdtVh 

Q 

“b V x'^h^ “b OhEx'^h * ^ x^t'^h) dxdt^ 


lh{vh) = J f[vh + OhdtVh] dxdt, 

Q 


(3.12) 

(3.13) 


and 0 is a positive constant which will be determined later. 

Next we will show that the discrete bilinear form ()3.12p is Vo/j—coercive 
with respect to the discrete norm 


h — 




2 

L2(St) 


1/2 


(3.14) 


Remark 1. We note that (j3.14l) is a norm on Vofe- Indeed, if ||u/j||/i = 0 for 
some Vh G Vo^, then we have ExVh = 0 and dtVh = 0 in Q and also Vh = 0 
on Et- This means that the function Vh is constant in the space-time 
cylinder Q. Since Vh = 0 on the boundary Et and all functions of the 
discrete space Vq/i are continuous, we conclude that = 0 in the whole 
space-time domain Q. The other norm axioms are obviously fulfilled. 
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Lemma 1. The discrete bilinear form ah{-,-) ■ Voh x Voh 1^, defined 
by (I3.12p . is Voh—coercive wrt the norm || • ||/i, i.e. 

ah{vh,Vh) > TchhWh '^Vh^Vok (3.15) 


with yLc = 1. 


dxdt 


Proof. Using Gauss’ theorem and the fact that nt = 0 on U, we get 

+ V... . . VMV,) dxdt 

Q 

= j ^dtvl + eh{dtVhf + \VxVh\^ + ^dt\Va,Vj 

Q 

1 /* Oh f 

=- J vlmds + eh\\dtVh\\l^(^Q) + \\vxVh\\l^(^Q) + Y J 

dQ dQ 

=l {hh\\l,^sr) - ll^/*llL(i;o)) +^^ll^i^'^lli2(Q) + 


VxVhl'^ntds 


for all Vh G Voh- Since Vh{x,0) = 0 and, therefore, Va;n/i(x,0) = 0 for all 
X £ 17, i.e., Vh = 0 and V^x/i = 0 on Uq, we immediately arrive at the 
estimate 

1 Oh 

C'hiVh, Vh) =2 ll^^lli2(X'T) ll^3;l^ft|li2(Q) + “^ll^3;''^/l|li2(I7T) 

= \\vh\\l Vx e Voh, 


which proves our lemma. □ 

Lemma [1] immediately yields that the solution Uh £ Voh of (13.111) is 
unique. Indeed, let us assume that there is another solution Uh £ Vq/i such 
that ah{uh,Vh) = lh{vh) for all Vh £ V^h- Then, taking the difference, we 
get ah{uh - Uh,Vh) = 0 for all Vh £ Voh- Choosing Vh = Uh - Uh ^ Voh 
and applying Lemma [U we immediately obtain the estimate 0 < \\uh— 
UhWh — c^h{uh — Uh,Uh — Uh) = 0, i.e., Uh — Uh = 0. Thus, the solution is 
unique. Since the discrete variational problem (j3.1ip is posed in the finite 
dimensional space Voh, the uniqueness yields existence of the solution 
Uh £ Voh of p.lip . 
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Since we already choose the standard B-spline resp. NURBS basis of 
Voh in (|3.9I) as Vq/i = span{(/j/i^i}jgx, we look for the solution Uh G Vbh of 
(I3.11h in the form of 

Uh{x, t) = Uh{xi, . . . , Xd, Xd+l) = ^ Ui(ph,i{x, t) 

i£X 

where Uf^ := [uiji^x G is the unknown vector of control points 

defined by the solution of the system of algebraic equations 

KhUh=[_h' 

with the Nh x system matrix Kh = [RTij = ah{iphj, ^h,i)]i,jei and the 
right-hand side = [/i = lh{'^h,i)]i€X G that can be generated as 
in a usual IgA code for elliptic problems. It is again clear from Lemma [1] 
that the stiffness matrix is regular. 

In the next section, we derive an a-priori discretization error estimate 
wrt the II • ||/i norm. 

4 Discretization Error Analysis 

At first, we will present the main ingredients, which are necessary for the 
derivation of our a priori discretization error estimate, in form of Lemmas. 


Lemma 2. Let K = <L>{K) G ICh and K G ICh- Then the scaled trace 
inequality 


\\v\\L2{dK) < Ctrhj^^"^ (||i^IIl2(A) + hK\v\m(K)) (4-1) 

holds for all v G H^{K), where hx denotes the mesh size of K in the phys¬ 
ical domain, and Ctr is a positive constant that only depends on the shape 
regularity of the mapping <P and the constant C in the trace inequality 

II/IIl2(9(o,i)‘*+i) — C'||/||j7i((o,i)d+i) (4.2) 

that is valid for all f G 

Proof. The proof follows from (|4.2p by mapping and scaling similar to 
the proof of m Theorem 3.2], where (14.11) is proved for d = 2 in such a 
way that the structure of the constants Ctr is explicitly given. □ 
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Lemma 3. Let K he an arbitrary mesh element from JCh- Then the in¬ 
verse inequalities 

\\^Vh\\L2{K) ^ Cinv,l hK\\Vh\\L2lK) (4-3) 

and 

\\'>Jh\\L2{dK) < Cinvfihp}^^\\Vh\\L2{K) ( 4 - 4 ) 

hold for all Vh G Vh, where Cinv,i and Cinv,o are positive constants, which 
are independent of h^ and K G IC^ . 


Proof. The inverse inequality (14.3p is a special case of the inverse inequal¬ 
ities given in [71 Theorem 4.1] and [71 Theorem 4.2]. The proof of (14.4h 
can be found in |12l Theorem 4.1]. □ 

We note that estimate (14.31) immediately implies the inverse inequality 

— Cinv,lhj,^ llL2(i<') G Vh, 

because dt is a part of V = {Vx,dt)~^. Below we need the inverse inequality 

\\dtdxi'ah\\L2{K) — Cinv,lhj^ II^Xi'l’h||L2(il') (4-5) 


that is obviously valid for all Vh G 14 and 1 = 1,..., d as well. 

To prove the a priori error estimate, we need to show the uniform 
boundedness of the discrete bilinear form ah{-,-) on Voh,* x Voh, where 
the space V^h,* = H -|- Vq/i is equipped with the norm 


I 11 h,* — 


H'i + {0h) ^WvhW'i^^Q) 


1/2 


(4.6) 


Lemma 4. The discrete bilinear form ah{-,-), defined by (I3.12p . is uni¬ 
formly bounded on Voh,* x Voh, he., there exists a positive constant pLi, that 
does not depend on h such that 

\ah{u,Vh)\ < fab\\u\\h,*\\vh\\h, Vu G Voh ,*,^vh G Vqh- (4.7) 

Proof. Let us estimate the bilinear form (j3.12p 

...(u,+ 0ha.ua,V, + v.u. V... + dun, 

Q 

term by term. 
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For the first term, since Vq/i C we can proceed with an integra¬ 
tion by parts wrt t, and then estimate the resulting terms by Cauchy- 
Schwarz inequality as follows: 


Q Q St 


< 


“ 0h\\dtVh\\l^^Q) 


+ 


Ml2{St) 


2 

L2(St) 


The second and third term of (j3.12D can easily be bounded by means 
of Cauchy’s inequality as follows: 


eh J dtudtVhdxdt < dh\\dtu\\l^f^Q^ " 0h\\dtVh\\'i^{Q) 

Q 


and 


/ 

Q 


VxU • Va;U/i dxdt < 




I'^xVhl 


L2(Q) 


The final term in the bilinear form is bounded from above by applying the 
Cauchy-Schwarz inequality, inverse inequality (14.51) . and inequality (13.71) 
to obtain 


6h J VxUh ■ xdtVh dxdt < 
Q 




{Ohf\\dtV:rVh\\l^^Q) 


< 


< 


IVtuI 


IVtuI 


IVtuI 


L2{Q) 


L2{Q) 


L2{Q) 


i=l 


i0h)'^C?nv,lYl ^KW^^i'^hWl^iK) 

i=i KeKh 


^ 2/^2 / 32 || y 7 ||2 
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Combining the terms from above and using Cauchy’s inequality, we get 


ah{uh,Vh) < 


+ 

X 


+ 


{dh) ^lklli2(Q) + ||«|li2(i;T) 

1 

+ ||V3;u|||2(Q) + ll^xit|li2(Q) 

1 


< 

X 

< 


2|| Va;U|| j;^2((3) + Oh\\dtu\\\^(^Q-^ + {9h) ^||M||i2(Q) + II'*^IIl2(i:'t)J 

ll'^ll h,* II II h ; 


where Hb = (2max{(l + I2, 1})^/^ □ 

Next we recall some approximation properties of our B-spline resp. 
NURBS space that follow from the approximation results proved in [71 
Section 3] and [171 Section 3]. Indeed, there exists projective operators 
nh '■ L 2 {Q) —)• Vh that deliver the corresponding asymptotically optimal 
approximation results. 


Lemma 5. Let I and s be integers with 0<l<s<p+l, and let 
V G H^{Q) n Hq'q{Q). Then there exist a projective operator 11h from 
Hq’q{Q) n H^{Q) to Voh and a positive generic constant Cg such that 

\n — < Cgh ^ ^II^IIh®(q)’ (4-8) 

KeKh 

where h again denotes the mesh-size parameter in the physical domain, p 
denotes the underlying polynomial degree of the B-spline resp. NURBS, 
and the generic constant Cg only depends onl,s andp, the shape regular¬ 
ity of the physical space-time domain Q described by the mapping <L> and, 
in particular, the gradient VL> of the mapping T>, but not on h and v. 

Proof. The proof follows the proof of the approximation results presented 
in Subsections 3.3 and 3.4 of [7|, see also [471 Proposition 3.1]. □ 
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If TlhV belongs to Vbh H H\Q) (the multiplicity of the inner knots is not 


larger than p -\- 1 — /), then estimate (|4.8I) immediately yields the global 
estimate 




(4.9) 


The basic approximation estimates (j4.8li and (I4.9jl yield estimates of 
the approximation error v — 11 hV wrt to the || • \\L 2 {dQ) - norm as well 
as wrt to the discrete norms || • ||/i and || • ||h,* which we later need in 
estimation of the discretization error u — Uh- 

Lemma 6 . Let s be a positive integer with 1 < s < p -|- 1, and let v G 
H^{Q) n i^Q q((5). Then there exist a projection 11 h from Hq’q{Q) to V^h 
and generic positive constants Ci, C 2 and C 3 such that 


V - nhv\\L^^9Q) <Cih^ ^^"^WvWhs^q), 

llu - IlhvWh < 

\\v — UhvWh,* < C^h^ 


(4.10) 

(4.11) 

(4.12) 


where h is the mesh-size parameter in the physical domain, p denotes 
the underlying polynomial degree of the B-spline resp. NURBS, and the 
generic constants Ci, C 2 and C 3 only depends on I, s and p, the shape 
regularity of the physical space-time domain Q described by the mapping 
<P and, in particular, the gradient of the mapping <L, but not on h and 

V. 

Proof. By applying inequality m, the quasi uniformity assumption 
(|3.7p and Lemma O the proof of ()4.10p is obtained as follows: 


V - UhvWl^^QQ-^ — ^ 11 ^ “ 


KeKh 

dKnaQ^% 


Y1 ^hV\\l^^K)+hK\v- 


KGiCh 



ag/c^ 



< 2Cl (Cuh-^Csh^^ + II 

< 2ClCs{Cu + l)/i2^-iu|||,.(Q) = Cik 



2s-l 
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The definition (|3.14p of the norm || • ||/i, the approximation error estimate 
(I4.8P for / = 1, and the estimate (I4.10p just proved yield 

11^ - nhv\\l =||v^(u - nhv)\\l^^Q) + eh\\dtiv - nhv)\\l^^Q) 


+ 2 !!^ 


iH-iQ) 


={Cs + Cs9h + = C 2 h^^^-^^\\v\\l 


which proves the second estimate of Lemma [6l Now let us prove the last 
estimate. Using the definition of the norm (14.6p . the just proven estimate 
P4.11I) . and the approximation error estimate (14.8p for / = 0, we get 


\v - = ||u - Uhvll,^ + {9h) ||u - HhvW 

(C 2 h^i^-^'> + {9h)-^Csh^^') \\v\\hiQ) 

|2 

iH-iQ) 


< 




= {C 2 + 9-^hCs) h^^^-^^\\v\\ls(o) = C3h2(^-^)||u|'2 

which completes the proof of the lemma. □ 

Lemma 7. If the solution u G Hq'^{Q) of the variational problem p2.1h 
belongs to then it satisfies the consistency identity 

ah{u,Vh) = lh{vh) V Ufe G Voh. (4.13) 

Proof. Since u G Hq’^{Q) n integration by parts in (|2.ip wrt t 

and X gives the variational identity 


/ fvdxdt+ / uqv ds = {dfU — Au)v dxdt — / ntuv ds + / nx-'^xuvds 
Q Eo Q dQ dQ 

= J {dtu — Au)v dxdt — J uv ds + j uv ds + J Ux ■ Vxuv ds 

(4.14) 


Eo 


dQ 


Q 

= J{dtu — Au)v dxdt + uv ds Vu G ^^^^(Q). 


Q 


Eo 


If we take a test function v G C^{Q) C then it follows that 


0,0' 


Q 


j(8,u-Au)vd.it. 


Q 
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Since C^{Q) is dense in L 2 {Q), we have 


dtU-Au = f in L 2 {Q). 


From (j4.14p . we now conclude that 


J uqv ds 

So 



uv ds 


for all V G i.e., u = uq on Uq. Since u G Hq (Q), the Dirichlet 

trace of tt on Z" is zero. This gives us the strong form of our model problem 

dni). 

Now, multiplying dtu — Au = / by a test function Vh + OhdtVh, where 
Vh £ Voh, integrating over Q, we get 

j{dtu - Au){vh + OhdtVh) dxdt = J f{vh + OhdtVh) dxdt 'ivh G Voh- 

Q Q 

Integrating by parts wrt x gives 

ISMVU + «hB,vu)i.At + + 6 h8,v,) 

Q 

- / n. ■ = / f(v, + dkBv,) d.it. 

dQ 


Q 


Using the fact that Uj, = 0 on both Uq and Ut, and Vh G Vq/o we have 

I a.u(v, + 0ka.v,) d.it + v,u . + 6ha.v,) dMt 

Q 

= j fivh + OhdtVh) dxdt, 


Q 


which completes the proof. □ 

Now we are in the position to prove the main result for this section, 
namely the a priori discretization error estimate in the discrete norm \\-\\h- 
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Theorem 1. Let u G Hq'^{Q) n H^{Q) with s >2 be the exact solution 
of our model problem and let Uh G Voh be the solution to the IgA 

scheme (j3.11j) . Then the discretization error estimate 

\\u - UhWh < (4.15) 

holds, where C is a generic positive constant, t = min{s,p+ 1}, and p 
denotes the underlying polynomial degree of the B-splines or NURBS. 

Proof. Subtracting the IgA scheme 

Vh') — Ihiphf G VqH 

from the consistence identity 

ah{u,Vh) = lh{vh), Vu/i G Voh, 

we obtain the so-called Galerkin orthogonality 

ah{u - Uh,Vh) = 0, '^VhGVoh, (4.16) 

that is crucial for the discretization error estimate. 

Using now the triangle inequality, we can estimate the discretization 
error u — Uh as follows 

\\u - Uh\\h < \\u - Nhullh + WHhU - Uh\\h- (4.17) 

The first term is nothing but the quasi-interpolation error that can easily 
be estimated by means of Lemma 0 The estimation of the second term on 
the right-hand side of ()4.17l) is more involved. Using the fact that UfiU — 
Uh G Voh, the Voh-6hipticity of the bilinear form ah{-,-) as was shown in 
Lemma [U the Galerkin orthogonality (I4.16P , and the boundedness of the 
discrete bilinear form Lemma [H we can derive the following estimates 

PcWNhU - UhWl < ahiLIhU - Uh, TlhU - Uh) 

= ahiLIhU - u, IlhU - Uh) 

< fJ.b\\nhU - u\\h,*\\nhu - Uh\\h- 


Hence, we have 


Wiihu - uhWh < {hb/hc)\\nhu - u\\h,*. 


( 4 . 18 ) 
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Inserting (j4.18jl into the triangle inequality (j4.17p and using the estimates 
(14.121) and (14.111) from Lemma [6l we have 

\\u - Uh\\h < \\u - nhu\\h + WHhU - Uh\\h 

< \\u - UhuWh + {fib/f^c)\\nhu - u\\h,* 

< Cs{l + fib/ ^||t;||jys(Q), 

which proves the discretization error estimate 14.151 with C = (^^(l -|- 

flb/fic)- □ 

5 Moving Spatial Computational Domain 

In this section, we will formulate the space-time scheme for a moving 
spatial domain Q := {{x,t) € S (0,T)}, where I7(t) C 

for d = 1, 2,3. 



Fig. 2. Moving spatial domain t G [0, T]. 


Let us now assume that the underlying polynomial degree p > 2 and 
the multiplicity m of all internal knots is less than or equal p — 1, i.e., 
Voh is always a subset of C^{Q). Thus, the IgA (B-spline or NURBS) 
space Voh is a subspace of Hq’q{Q) n H^{Q). Similar to Subsection 13.21 
we derive the IgA scheme by multiplying our parabolic PDE (II.ip by a 
test function of the form Vh + OhdfVh with an arbitrary Vh G Voh) where 
d is a positive constant which will be defined later, and integrating over 
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the whole space-time cylinder Q, giving the identity 

J f{vh + OhdtVh) dxdt = j{dtUhVh + OhdtUhdtVh (5.1) 


Q 


Q 


- Auh{vh + OhdtVh)) dxdt 


that is valid for all Vh G Vbh- 

Integration by parts with respect to x in the last term of the bilinear 
form on the right-hand side of (|5.1I) gives 

J fivh + OhdtVh) dxdt = j{dtUhVh + OhdtUhdtVh + VxUh ■ VxVh 

Q Q 

+ OhVxUh ■ ^xdtVh) dxdt - j rix ■ {vh + OhdtVh) ds. 

dQ 

We now integrate the fourth term on the right-hand side by parts wrt 
time t and obtain 

jidm + 6h8,uAvH + V.Uk . - 6h8,VM ■ V.%) 

Q 

+ / . V... + eh8,v,) *. 


Q 


dQ 


dQ 


Using the facts that Vh G Voh and Hx is zero on Sq U Ut, we can 
continue to write 

= jia.u,v, + 6h8.uAv, + - 6h8,V,u, ■ i.At 

Q 

Considering the terms on the boundary U, we note the following 


Oh 


i: 
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since [ntVxVh — nxdtVh] is the tangential derivative of Vh and t;/! = 0 on 
E. Finally, we arrive at our discrete scheme: find Uh G Vbh such that 

bhiuh,Vh) = lh{vh) yvh G Voh, (5.2) 


where the bilinear and linear forms are given as follows: 


bh{uh,Vh) ■= j (dtUhVh + OhdtUhdtVh + VxUh ■ ^xVh 

Q 

- ewv.w ■ dxdt + ehj vm ■ v.% (s.s) 




lh{vh) '■= J f + OhdtVh] dxdt. 
Q 


Remark 2. Note that, for Uh,Vh G Vbh with p >2, we obtain an alterna¬ 
tive representation for the discrete bilinear form 

h{u,,v,) = + ekS,u,Av, + 

Q 

+ . W8,v,)dxdt + 0hjn,y.uu . V... * 

i: 

by integration by parts. For non-moving domains, we observe that re* = 0 
on 17, and we therefore obtain 

bhiuh, Vh) = ah{uh,Vh) for all Uh, Vh G Voh- 

Hence, in the case of a non-moving domain, the discrete variational for¬ 
mulation (j3.1ip is equivalent to the discrete problem (15.2p . 

In the next lemma, we show that the discrete bilinear form for the 
moving spatial domain is Vq/i— coercive wrt to the energy norm 

1 

v\\h~^ dhWS/ 



Lemma 8. Let 9 > 0 be sufficiently small, in particular, we assume 

9 < i2Cinv,lCu)~\ 
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where Cinv,i is the constant from Lemma 0 and the constant Cu > 0 is 
given by the quasi uniform assumption Then the discrete bilinear 

form bh{-,-) ■ Voh X Voh M, defined by (15.31) . is Voh—coercive wrt the 
norm || • \\h,m, i-c. 

bh{vh,Vh) > hc\\vh\\h,m^ Vu G Vb/i (5.4) 

with p,(. = 1/2. 

Proof. Let Vh G Vok- By means of the definition of the discrete bilinear 
term bh{-, ■) from (|5.3p and using Gauss’ theorem, we get 

Q 


Q 


dtvl + eh{dtVhf + \S/xVh? - OhdtVxVh ■ ^xVh 


dxdt 


+ 0/1 y \VxVh\ds 


=\ j clntds + eh\\dtVh\\l^^Q) + \\VxVh\\l^^Q) 

dQ 

Q 

for all V G V^h- Since r;/i = 0 on 27 U i7o, we further obtain 

bh{vh,Vh) =^\\vh\\l^(ET) + +^^ll^=^^^lli2(i:T) 

Q 

>\\vh\\l,m - GhWdtVxVh\\L 2 {Q)\\^xVh\\L 2 {Q) ^r all Vh G Vqh- 

By using the inverse inequality (|4.5I) and the assumption (13.7p we further 
obtain similar to the proof of Lemma |4] the estimate 


(I 


L2{Q) ~ X] xVh\\~L.^(^K) - 2-^ xVhWf.^iK) 

K&Kh K&Kh 


|2 (-i 2 /i- 2 ||y 7 Il 2 


— ^inv,l^u^ ^ xVhW'i.^iK) ~ [Cinv,lCuh ^ || || ^^(Q)] • 

KeiCh 




Space-Time IgA of Parabolic Evolution Equations 


23 


Hence, we have for 6 < {2Cinv,iCu) ^ the estimate for the bilinear form 

bhivh,Vh) > \\Vh\\l,m-^HdtVr,Vh\\L 2 {Q)\\'^xVh\\L 2 (Q) 

> (1 - ehCinv,lCuh-^)\\vh\\l,^ > ^\\vh\\l,m- 

□ 


To show the boundedness for the bilinear form bh{-,-), we need the fol¬ 
lowing additional norm 




\v\\h,m + idh) ^lln 


|2 

'L 2 {Q) 


+ i9^h^)\\dtV,v 


|2 

\L 2 (Q) 


Lemma 9. The discrete bilinear form bh{-,-) ■ ^ defined 

by (15. is uniformly bounded on Vq/i,* x *-e., there exists a positive 
constant yn, which does not depend on h such that 


I (tt, n/j) I ^ llrtIIIl't’/i II/i,m (^■^) 

for all u G Voh,* and all Vh G Voh- 
Proof. We estimate the discrete bilinear form ()5.,‘lj) 


bh{u, Vh) = j {dtuvh + OhdtudtVh + ^xU ■ VxVh - OhdtVxU ■ VxVh) dxdt 


Q 


term by term. For the first, second and third term, we can proceed as in 
the non-moving case, and we obtain the estimates 


J dtuvhdxdt < (eh) 2 ||m||l2(q)(6»/i) 2 ||9tn/,||i2(Q) 

Q 


+ \W\\L2{i;T)\\'^h\\L2{ST)^ 

eh J dtudtVhdxdt < {eh)^\\dtu\\L^(^Q){0h)^\\dtVh\\L2{Q}, 

Q 

j V^^u-V^Vhdxdt < \\yxU\\L2{Q)\\^xVh\\L2(Q)- 

Q 
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For the two remaining terms, we also use the Cauchy-Schwarz inequality, 
and we obtain 


9h j dtVxU-VxVhdxdt < eh\\dtVxu\\L^(^Q)\\Va^Vh\\L2{Q), 

Q 

Oh J -VxVhds < {9h)^\\Vxu\\L^(^ST)i^h)^\\'^xVh\\L2{i:T)- 

X'x’ 

Combining these estimates, we conclude the statement of this Lemma 
with 


bh{u,Vh)< {9h) ^\\u\\l^^Q) + + 9h\\dtu\\l^(^Q^ 

+ W'^ xu\\'i^(^Q-j + 9'^h‘^\\dtVxU\\‘j^^(^Q^ + 

X dh\\dtVh\\l^^Q) + \\vh\\l^(ET) + ^f^W^tVhWl.^Q) 

+ xVh\\%(Q) + xVh\\%(Q) + dh\\V xVh\\‘i^(^ST) 

— 2||rillllu/i II 


□ 

Lemma 10. Let s be a positive integer with 2 < s < p + 1, and let 
V G H^{Q). Then there exist a projection Ilh from Hq’q{Q) n H^{Q) 
to Voh, and generic constants Ci,C 2 > 0 such that the following error 
estimates hold 


||u - 77/,u||fe,m < ^lblk'>(0)! (5-6) 

lb - nhv\\h,m,* < C'2h^“bbll//'>(Q)> (5-7) 

where p denotes the underlaying polynomial degree of the B-spline resp. 
NURBS. 

Proof. By using Lemma [6] it remains to estimate the additional terms 
h^\\dtVx{v-nhv)\\l^(^Q^ and h||Va;(u - 77,,u)|||^(^^). 


v\\ho(Q)- 


With Lemma E] we obtain for i, j = 1,..., d + 1 
Wdxidx^iv - 77/,u)|b2(Q) < Csh^~^ 


(5.8) 
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By using the error estimate (|5.8p we obtain 

d 


h?\\dtyx{v - nhv)\\l^^Q^ = Wdtdxiiv - nhv)\\l^^Q^ 

i=\ 

i=\ 

VdCgh^ ^ll^llff»(Q) 

With Lemma [5] and using also (|5.8I) we further obtain 


h||Vx(u nhv)\\\^(^x:T) ~ ^ X!/ -^/i^)llL2(9Ani:T) 


KeiCh *=i 
dKnSTj^m 


< 


'> E E^J hx - nhV)\\L2{K) + hK\d:,,{v - il/jU)|' 


K&Kh i=l 
dKnEr^ili 


a 

< 2 CuCl Y. E - nhv)\\YK) + h^\d-Ay - nhv)\l.(K)) 


K&Kh i=l 
dKnST^tl) 


< 2 CuCl 


= 2 CuCt 


|Va;(u - 77/,u)||i^(Q) + h^Y 


(Q) 


2=1 

d d-\-l 


|v.(«-n»q|lL,g, + ft^EE \dxidxjiv - nhv)\jji 


iQ) 


< 


^/2{l + d{d + l))CuCtrCsh^-^ llu 


i=l j=l 
1 ^ 

\h‘^(Q) 


which completes the proof. 


□ 


Lemma 11. Let p > 2. If the solution u G Hq’^{Q) of the variational 
problem m belongs to H'^{Q), then it satisfies the consistency identity 


bhiu, Vh) = lh{vh) V u/i € Voh- 


Proof. With the same arguments as in Lemma [71 we obtain that 

dtu — Au = f in L 2 (Q) and u = uq in ^ 2 (^ 0 ). (5-9) 
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We now multiply the differential equation of (j5.9|) with a test function 
Vh + OhdtVh for Vh G Vq/i and integrate over the space-time domain Q. 
Since u € H^{Q) and p> 2, we can apply all the derivations as we did at 
the beginning of this section to obtain the statement of this Lemma. □ 

Theorem 2. Let p > 2 and 6 be sufficiently small, see Lemma 0 Fur¬ 
thermore, let u G Hq'^{Q) n H^{Q) with s > 2 be the exact solution of 
the model problem (EH), and let Uh G Voh be the solution to discrete 
variational problem dSH- Then the discretization error estimate 

Ch^ ^\\u\\Ht{Q), (5.10) 

holds, where t = min{s,p-|- 1}, and p denotes the underlying polynomial 
degree of the NURBS. 

Proof. The stated error estimate follows in the exact way as in Theorem 
mby using Lemma IHHITl □ 

Remark 3. For the case p = 1, Lemma [IT] is not valid in general. In 
this case, we have to take additional consistency errors into account by 
estimating them properly. It turns out, cf. Subsection 16.21 that also, for 
the case p = 1, we obtain the full convergence rates wrt the energy norm 
II • ||h,m for smooth solutions. 

6 Numerical Results 

The numerical results presented below have been performed in G+SMC0 
[23]. We used the sparse direct solver SuperLU to solve the resulting lin¬ 
ear system ()3.16p of IgA equations. We present numerical results for both 
fixed fSubsection 16. 1 p and moving fSubsection l6.2p spatial computational 
domains in one and two dimensions in space. In Subsection [631 we present 
numerical results which demonstrate the efficiency of a standard paral¬ 
lel AMG preconditioned GMRES solver on massively parallel computers 
with several thousands of cores. We mention that 0 was set to 0.1 in all 
our numerical experiments including the examples with moving spatial 
domains. 

6.1 Fixed Spatial Computational Domain 


^ Geometry plus Simulation Modules, http://www.gs.jku.at/gismo 
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6.1.1 Fixed one-dimensional spatial computational domain. We 

consider the one-dimensional spatial domain 1? = (0,1) and the time in¬ 
terval (0, T) with r = 1, i.e., we have the space-time cylinder Q = (0,1)^, 
that can geometrically represented by the knot vectors Si = {0,0,1,1} 
and S 2 = {0,0,1,1} in the IgA context. We solve our parabolic boundary 
value problem (11.11) . and choose the data such that the solution is given 
by u{xi,t) = u{x, t) = sin(7rx) sin(7rt), i.e. f{x, t) = dtu{x, t) — Au{x, t) = 
TT sin(7rx)(cos(7rt) -|- 7rsin(7rt)) in Q, uq = 0 on 17, and u obviously van¬ 
ishes on S. Thus, the compatibility condition between boundary and 
initial conditions holds. The convergence behavior of the space-time IgA 
scheme with respect to the discrete norm || • ||/i is shown in Tables [1] and [2] 
by a series of /i-refinement and by using B-splines of polynomial degrees 
p = 1,2,3,4. After some saturation, we observe the optimal convergence 
rate 0{hP) as theoretically predicted by Theorem [T] for smooth solutions. 
Moreover, Tables [3] and 0] show the L 2 errors and the corresponding rates 
for the same setting. We see that the L 2 rates are asymptotically optimal 
as well, i.e. the L 2 -error behaves like 



0.0 0.2 0.4 0.6 0.8 1.0 

x-Axis 


Fig. 3. Solution contours in the space-time cylinder Q for Example 16.1.11 
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Table 1. E 

and p = 2. 


Table 2. Errors and rates wrt 
and p = A. 


p = 1 

p = 2 

Dofs 

W - Uh\\h 

Rate 

Dofs 

||a - Uh\\h 

Rate 

4 

1.6782e-H00 

0 

9 

2.10399e-01 

0 

9 

7.28214e-01 

1.20 

16 

2.03729e-01 

0.04 

25 

3.61278e-01 

1.01 

36 

3.98228e-02 

2.35 

81 

1.79489e-01 

1.01 

100 

9.29436e-03 

2.10 

289 

8.94084e-01 

1.01 

324 

2.27848e-03 

2.02 

1089 

4.46132e-02 

1.00 

1156 

5.66197e-04 

2.01 

4225 

2.22829e-02 

1.00 

4356 

1.41258e-04 

2.00 

16641 

1.11354e-02 

1.00 

16900 

3.52865e-05 

2.00 

and rates wrt 

||/l 

for Example 1 6.1 

.l|ar 

p = 3 

P = 4 

Dofs 

M — Uh\\h 

Rate 

Dofs 

U - Uh\\h 

Rate 

16 

2.10106e-01 

0 

25 

6.70957e-03 

0 

25 

2.73234e-02 

2.92 

36 

6.49630e-03 

0.05 

49 

5.08124e-03 

2.43 

64 

5.56376e-04 

3.55 

121 

5.73528e-04 

3.15 

144 

3.27142e-05 

4.10 

361 

6.93807e-05 

3.05 

400 

2.05481e-06 

4.00 

1225 

8.58843e-06 

3.01 

1296 

1.30057e-07 

4.00 

4489 

1.07029e-06 

3.00 

4624 

8.20252e-09 

4.00 

17161 

1.33647e-07 

3.00 

17424 

5.15378e-10 

4.00 


for Example 16.1.11 and degrees p = 3 


P = 1 

p = 2 

Dofs 

||m - Uh\\L2(Q) 

Rate 

Dofs 

\\u-Uh\\L2(Q) 

Rate 

4 

5.000e-01 

0 

9 

2.69186e-02 

0 

9 

1.21333e-01 

2.04 

16 

2.66817e-02 

0.03 

25 

3.03720e-02 

2.00 

36 

2.30767e-03 

3.53 

81 

7.47639e-03 

2.02 

100 

2.60187e-04 

3.15 

289 

1.85552e-03 

2.01 

324 

3.16609e-05 

3.04 

1089 

4.62271e-04 

2.00 

1156 

3.92785e-06 

3.01 

4225 

1.15377e-04 

2.00 

4356 

4.89712e-07 

3.01 

16641 

2.88212e-05 

2.00 

16900 

6.11484e-08 

3.00 


Table 3. L 2 errors and rates for Example 16.1.11 nsing degrees p = 1 and 

p = 2. 


6.1.2 Fixed two-dimensional spatial computational domain. As 

a second example we consider the two-dimensional spatial domain 17 = 
(0,1)^ and the time interval (0, T) with T = 1, i.e., we have the space-time 
cylinder Q = (0,1)^, that can geometrically be represented by the knot 
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p = 3 

p = 4 

Dofs 

W - Uh\\L2{Q) 

Rate 

Dofs 

\\u-Uh\\L2(Q) 

Rate 

16 

2.65068e-02 

0 

25 

5.49027e-04 

0 

25 

2.32514e-03 

3.16 

36 

5.43598e-04 

0.01 

49 

3.10354e-04 

2.84 

64 

3.90240e-05 

3.80 

121 

1.63884e-05 

4.25 

144 

1.02403e-06 

5.25 

361 

9.72892e-07 

4.10 

400 

3.03989e-08 

5.07 

1225 

5.99946e-08 

4.02 

1296 

9.38131e-10 

5.02 

4489 

3.73705e-09 

4.01 

4624 

2.92231e-ll 

5.00 

17161 

2.33370e-10 

4.00 

17424 

9.15973e-13 

5.00 


Table 4. L 2 errors and rates for Example 16.1.11 using degrees p = 3 and 
p = 4. 


vectors Si = {0,0,1,1}, S2 = {0,0,1,1} and S’3 = {0,0,1,1} in the con¬ 
text of IgA. We solve our model problem (ll.ljl . and again choose the data 
such that the solution is given by u{x,t) = sin(7ra:i) sin(7ra:2) sin(7rt), i.e. 
f{x, t) = dtu{x, t) — Au{x, t) = 7rsin(7rxi) sin(7ra:2)(cos(7rt)-|-27r sin(7rt)) in 
Q, uq = 0 on 17, and u obviously vanishes on S. Thus, the compatibility 
condition between boundary and initial conditions holds. In Figure 01 we 
present the solution contours of the problem in where we have sliced 
the domain along the f-axis. The convergence behavior of the space-time 
IgA scheme with respect to the discrete norm || • ||/i is shown in Tables [5] and 
[6] by a series of h-rehnement and by using B-splines of polynomial degrees 
p = 1,2,3,4. After some saturation, we observe the optimal convergence 
rate 0{hP) as theoretically predicted by Theorem [1] for smooth solutions. 
Moreover, Tables [3 and [8] show the L 2 errors and the corresponding rates 
for the same setting. We see that the L 2 rates are asymptotically optimal 
as well, i.e. the T 2 -error behaves like 0{hP^^). 


Table 5. Errors and rates wrt 


p = 1 

P = 2 

Dofs 

M - Uh\\h 

Rate 

Dofs 

||u - Uh\\h 

Rate 

8 

1.63740 

0 

27 

2.15440e-01 

0 

27 

7.39981e-01 

1.15 

64 

2.11247e-01 

0.03 

125 

3.60495e-01 

1.04 

216 

3.98871e-02 

2.40 

729 

1.79065e-01 

1.01 

1000 

9.27926e-03 

2.10 

4913 

8.92787e-02 

1.00 

5832 

2.27556e-03 

2.03 

35937 

4.45779e-02 

1.00 

39304 

5.65772e-04 

2.01 


for Example 16.1.21 and degrees p = 1 


and p = 2. 
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Fig. 4. Solution contours in the space-time cylinder Q at t = 0.5 for 
Example 16.1.21 


Table 6. Errors and rates wrt 
and p = 4. 


p = 3 

p = 4 

Dofs 

M — Uh\\h 

Rate 

Dofs 

U - Uh\\h 

Rate 

64 

2.16883e-01 

0 

125 

6.67416e-03 

0 

125 

2.75120e-02 

2.97 

216 

7.69213e-03 

0.03 

343 

5.09465e-03 

2.43 

512 

3.55820e-03 

3.55 

1331 

5.72742e-04 

3.15 

1728 

3.26623e-05 

4.10 

6859 

6.92964e-05 

3.04 

8000 

2.05215e-06 

4.00 

42875 

8.58214e-06 

3.01 

46656 

1.37611e-07 

4.00 


for Example [67 


2] and degrees p = 3 


6.2 Moving Spatial Computational Domain 

6.2.1 A simple one-dimensional moving spatial computational 
domain. Now we consider the one-dimensional moving computational 
domain f2{t) = {x = xi £ : a{t) < x < b{t)} with t = (0,T), 

where a{t) = —tj2, b[t) = 1 -|- t/2, and T = 1. The space-time cylinder 
Q = {(x, t) G : X G G (0, T)} C is obviously a fixed domain 
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p = 1 

p = 2 

Dofs 

M - Uh\\L2{Q) 

Rate 

Dofs 

\\u-Uh\\L2(Q) 

Rate 

8 

3.65528e-01 

0 

27 

2.39153e-02 

0 

27 

9.56396e-02 

1.93 

64 

2.37388e-02 

0.01 

125 

2.32679e-02 

2.03 

216 

1.99848e-03 

3.57 

729 

5.75358e-03 

2.01 

1000 

2.22710e-04 

3.17 

4913 

1.43171e-03 

2.01 

5832 

2.70486e-05 

3.04 

35937 

3.57195e-04 

2.00 

39304 

3.35780e-06 

3.00 


Table 7. L 2 errors and rates for Examp 

p = 2. 


e 16.1.21 using degree p = 1 and 


p = 3 

p = 4 

Dofs 

||m - Uh\\L2(Q) 

Rate 

Dofs 

I|w-w'iIU2{q) 

Rate 

64 

2.39325e-02 

0 

125 

4.75391e-04 

0 

125 

2.03108e-03 

3.56 

216 

4.73425e-04 

0.01 

343 

2.68174e-04 

2.92 

512 

3.34666e-05 

3.82 

1331 

1.41715e-05 

4.24 

1728 

8.74291e-07 

5.26 

6859 

8.42223e-07 

4.07 

8000 

2.60544e-08 

5.07 

42875 

5.19528e-08 

4.01 

46656 

8.07120e-10 

5.01 


e 16.1.21 using degree p = 3 and 


p 


= 4. 


in the space-time world M^. It can geometrically be represented by the 
knot vectors Si = {0,0,1,1} and S 2 = {0,0,1,1} and the control points 
Pi,i = (0,0), P 2,1 = (1,0), P 2,2 = (1.5,1.0), and Pi ,2 = (-0.5,1.0) in 
the context of IgA. We solve our model problem (II.Ill , and again choose 
the data such that the solution is given by u{x,t) = sin(7rx) sin(7rt), i.e. 
f{x,t) = dtu{x,t) — Au{x,t) = (tt sin(7rx))(cos(7rt) -|- 7rsin(7rt)) in Q, 
uo = 0 on 17, and u{x, t) = sin(7ra;) sin(7rt) on S. Thus, the compatibility 
condition between boundary and initial conditions holds. The space-time 
computational domain Q and the solution is drawn in Figure O The 
convergence behavior of the space-time IgA scheme with respect to the 
discrete norm || • \\h,m is shown in Tables M and 1101 by a series of h- 
refinement and by using B-splines of polynomial degrees p = 1,2, 3,4. 
After some saturation, we observe the optimal convergence rate 0{hP) 
for p > 2 as theoretically predicted by Theorem [2] for smooth solutions. 
Moreover, Tables [7] and [8] show the L 2 errors and the corresponding rates 
for the same setting. We see that the L 2 rates are asymptotically optimal 
for p>2 as well, i.e. they behave like For p = 1, we also observe 

the optimal rate in the discrete norm, cf. also Remark [3l whereas the 
L 2 -rate does not reach the optimal order 2. 
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Fig. 5. Solution contours in the space-time cylinder Q for Example 16. 2. 11 


p = 1 

p = 2 

Dofs 

||w 

Rate 

Dofs 

11'^^ 'ti/i II 

Rate 

4 

2.31256 

0 

9 

8.3047e-01 

0 

9 

8.7222e-01 

1.41 

16 

7.4075e-01 

0.17 

25 

5.9385e-01 

0.55 

36 

1.3091e-01 

2.50 

81 

3.0387e-01 

0.97 

100 

2.9977e-02 

2.13 

289 

1.5240e-01 

1.00 

324 

7.2872e-03 

2.04 

1089 

7.6181e-02 

1.00 

1156 

1.8068e-03 

2.01 

4225 

3.8070e-02 

1.00 

4356 

4.5056e-04 

2.00 

16641 

1.9028e-02 

1.00 

16900 

1.1255e-04 

2.00 

rs and rates wrt 

1 h,m 

or Examp 

e|6.:. 


16.2.11 and degn 


ees 


p = 1 and p = 2. 


6.2.2 Curvilinearly moving one-dimensional spatial computa¬ 
tional domain. We consider again an one-dimensional moving spatial 
computational domain of the form = {x = xi G : a{t) < x < 
b{t)}, t G (0,1), where now the movement is described by the functions 
a{t) = t{l — t)/2 and b{t) = 1 — t(l —1)/2 leading to the space-time cylin¬ 
der Q = {(x,t) G : X G G (0, T)} C with a curved surface 

area S, see Figure [6l The space-time cylinder Q can also be represented 
by the knot vectors Si = {0,0, 1 , 1} and S 2 = {0,0,0,1, 1 , 1} and the cor¬ 
responding control points Pi^ = (0,0), P 2 ,i = (1,0), P 2,2 = (0.75,0.5), 
P 2,3 = (1,1), Pi ,3 = (0,1) and Pi ^2 = (0.25,0.5) in the context of IgA, 
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Table 10. 

p = 3 and p = A. 


p = 3 

p = 4 


Dofs 


Rate 

Dofs 

II 

Rate 


16 

7.87827e-01 

0 

25 

9.58143e-02 

0 


25 

7.74683e-02 

3.35 

36 

7.29388e-02 

0.35 


49 

2.68888e-02 

1.53 

64 

7.25760e-03 

3.31 


121 

2.54304e-03 

3.40 

144 

3.25434e-04 

4.48 


361 

2.97843e-04 

3.09 

400 

1.83582e-05 

4.15 


1225 

3.69216e-05 

3.01 

1296 

1.11825e-06 

4.04 


4489 

4.62600e-06 

3.00 

4624 

6.95208e-08 

4.01 


17161 

5.79880e-07 

3.00 

17424 

4.34224e-09 

4.00 


irs and rates wrt 

|/i,m 

or Examp 

e|b.2 



]T] using degree 


p = 1 

p = 2 

Dofs 

l|w - Uh\\L2(Q) 

Rate 

Dofs 

l|w - 'Uh\\L2(Q) 

Rate 

4 

8.7462 

0 

9 

1.9450e-01 

0 

9 

1.9970e-01 

2.13 

16 

1.41565e-01 

0.46 

25 

7.3226e-02 

1.45 

36 

1.09044e-02 

3.70 

81 

1.8794e-02 

1.96 

100 

1.12045e-03 

3.28 

289 

5.0344e-03 

1.90 

324 

1.31179e-04 

3.09 

1089 

1.5998e-03 

1.65 

1156 

1.60564e-05 

3.03 

4225 

6.4746e-04 

1.31 

4356 

1.99240e-06 

3.01 

16641 

3.0437e-04 

1.10 

16900 

2.48337e-07 

3.00 


p = 2. 


p = 3 

p = 4 

Dofs 

||u - Uh\\L2{Q) 

Rate 

Dofs 

\\u-Uh\\L2(Q) 

Rate 

16 

1.69684e-01 

0 

25 

1.38383e-02 

0 

25 

1.02458e-02 

4.05 

36 

1.01827e-02 

0.44 

49 

2.34449e-03 

2.13 

64 

6.76183e-04 

3.91 

121 

9.90173e-05 

4.57 

144 

1.2387e-05 

5.77 

361 

5.56308e-06 

4.15 

400 

3.23506e-07 

5.26 

1225 

3.40489e-07 

4.03 

1296 

9.61797e-09 

5.07 

4489 

2.12529e-08 

4.00 

4624 

2.96841e-10 

5.02 

17161 

1.33061e-09 

4.00 

17424 

9.24929e-12 

5.00 


p = A. 


see also Figure [H We solve our model problem (11.11) . and again choose 
the data such that the solution is given by u{x,t) = sin(7rx) sin(7rt), i.e. 
f{x,t) = dtu{x,t) — Au{x,t) = (7rsin(7rx))(cos(7rt) -|-7rsin(7rt)) in Q, 
uq = 0 on 17, and u{x, t) = sin(7ra;) sin(7rt) on U. Thus, the compatibility 
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condition between boundary and initial conditions hold. The convergence 
behavior of the space-time IgA scheme with respect to the discrete norm 
II • \\h,m is shown in Tables [13] and [14] by a series of /i-refinement and by 
using B-splines of polynomial degrees p = 1, 2,3,4. After some saturation, 
we observe the optimal convergence rate 0{hP) for p > 2 as theoretically 
predicted by Theorem [2] for smooth solutions. Moreover, Tables [15] and 
M show the L 2 errors and the corresponding rates for the same setting. 
We see that the L 2 rates are asymptotically optimal for p > 2 as well, i.e. 
they behave like For p = 1, we also observe the optimal rate in 

the discrete norm, cf. also Remark [S] whereas the L 2 -rate does not reach 
the optimal order 2. 



Fig. 6. Solution contours for Example 16.2.21 (left) and the control points 
(right). 


6.2.3 Curvilinearly moving two-dimensional spatial computa¬ 
tional domain. We now consider the curvilinearly moving two-dimensional 
spatial computational domain f2{t) = {x = {xi,X 2 ) G : a{t) < xi < 
b{t),0 < X 2 < 1} in M^, where a{t) = t{l — t)/2, b{t) = 1 — t(l — t)/2, 
and t is running from 0 to T = 1, leading to the space-time cylinder 
Q = {(a;,t) G : x G Q{t),t G (0,T)} C that is fixed in the 
space-time world M^. In the context of IgA, Q can geometrically be 
represented by the knot vectors Si = {0,0,1,1}, S 2 = {0,0,1,1} and 
E’a = {0,0,0,1,1,1}, and the control points given in Figure [7] (right). 
We solve our model problem (II.ip . and again choose the data such that 
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Table 13. 

p = 1 and p = 2. 


p = 1 

p = 2 

Dofs 


Rate 

Dofs 


Rate 

4 

1.40608e-f00 

0 

9 

1.79646e-01 

0 

9 

6.28191e-01 

1.15 

16 

1.66950e-01 

0.11 

25 

3.02013e-01 

1.05 

36 

2.52637e-02 

2.72 

81 

1.47951e-01 

1.02 

100 

5.45976e-03 

2.21 

289 

7.33135e-02 

1.01 

324 

1.29401e-03 

2.08 

1089 

3.65079e-02 

1.01 

1156 

3.16214-04 

2.03 

4225 

1.82188e-02 

1.00 

4356 

7.82250e-05 

2.02 

16641 

9.10085e-03 

1.00 

16900 

1.94575e-05 

2.01 


\ h^m 


for Examp 


e 16.2.21 and degrees 


CO 

II 

a 

p = 4 

Dofs 

||w 

Rate 

Dofs 

ll^i 'U/i||^,7TT, 

Rate 

16 

2.25659e-01 

0 

25 

1.55287e-02 

0 

25 

1.95569e-02 

3.18 

36 

9.35760e-03 

0.73 

49 

3.37111e-02 

2.56 

64 

3.74209e-04 

4.64 

121 

3.37580e-03 

3.25 

144 

1.93936e-05 

4.27 

361 

3.95998e-04 

3.09 

400 

9.83348e-07 

4.30 

1225 

4.83941e-06 

3.03 

1296 

5.54882e-08 

4.14 

4489 

5.99023e-07 

3.01 

4624 

3.30365e-09 

4.07 

17161 

7.45335e-08 

3.01 

17424 

2.01770e-10 

4.03 

3rs and rates wrt 

1 h,m 

or Examp 

e|6.2 


16.2.21 using degree 


p = 3 and p = A. 


P = 1 

p = 2 

Dofs 

||m - Uh\\L2(Q) 

Rate 

Dofs 

\\u-Uh\\L2(Q) 

Rate 

4 

4.03831e-01 

0 

9 

2.93862e-02 

0 

9 

9.12333e-02 

2.14 

16 

3.21213e-02 

0.13 

25 

2.33973e-02 

1.96 

36 

2.74769e-03 

3.55 

81 

5.79121e-03 

2.01 

100 

3.10332e-04 

3.15 

289 

1.52473e-03 

1.92 

324 

3.77756e-05 

3.03 

1089 

4.68079e-04 

1.70 

1156 

4.68649e-06 

3.01 

4290 

1.82024e-04 

1.36 

4356 

5.84132e-07 

3.00 

16770 

8.37613e-05 

1.12 

16900 

7.29169e-08 

3.00 


p = 2. 


the solution is given by u{x,t) = sin(7rxi) sin(7r3:2) sin(7rt), i.e. f{x,t) = 
dtu{x,t) — Au{x,t) = (vr sin(7rxi) sin(7rx2))(cos(7rt) -|- 27rsin(7rt)) in Q, 
uq = 0 on i7, and u{x,t) = sin(7rxi) sin(7rx2) sin(7rt) on U. In Figure [7] 
(left), we present the solution contours of the problem in at t = 0.5. 
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p = 3 

p = 4 

Dots 

W - Uh\\L2(Q) 

Rate 

Dofs 

\\u-Uh\\L2(Q) 

Rate 

16 

3.22870e-02 

0 

25 

1.00347e-03 

0 

25 

2.61893e-03 

3.62 

36 

9.49495e-04 

0.10 

49 

4.05805e-04 

2.69 

64 

5.65065e-05 

4.10 

121 

1.99434e-05 

4.34 

144 

1.96066e-06 

4.85 

361 

1.15336e-06 

4.11 

400 

6.03483e-08 

5.02 

1225 

7.10583e-08 

4.02 

1296 

1.88123e-09 

5.00 

4489 

4.43879e-09 

4.00 

4624 

5.88213e-ll 

5.00 

17161 

2.77796e-10 

4.00 

17424 

1.83916e-12 

5.00 


p = A. 


The convergence behavior of the space-time IgA scheme with respect to 
the discrete norm || • \\h,m is shown in Tables and [18] by a series of 
/i-refinement and by using B-splines of polynomial degrees p = 1,2,3,4. 
After some saturation, we observe the optimal convergence rate 0{hP) 
for p > 2 as theoretically predicted by Theorem [2] for smooth solutions. 
Moreover, Tables M and [201 show the L 2 errors and the corresponding 
rates for the same setting. We see that the L 2 rates are asymptotically 
optimal for p > 2 as well, i.e. they behave like 0{h^^^). For p = 1, we 
also observe the optimal rate in the discrete norm, cf. also Remark (3] 
whereas the L 2 -rate does not reach the optimal order 2. 


p= 1 

p = 2 

Dofs 

11'^ 

Rate 

Dofs 

11'^ 

Rate 

8 

1.57742e-f00 

0 

27 

2.5829e-01 

0 

27 

6.84322e-01 

1.20 

64 

2.0891e-01 

0.31 

125 

3.25869e-01 

1.07 

216 

3.64503e-02 

2.52 

729 

1.60291e-01 

1.02 

1000 

7.88416e-03 

2.21 

4913 

7.96709e-02 

1.01 

5832 

1.85329e-03 

2.08 

35937 

3.97383e-02 

1.00 

39304 

4.50536e-04 

2.04 


Table 17. Errors and rates wrt 
p = 1 and p = 2. 


h.m 


:'or Example 16.2 


3] using degree 


6.3 Parallel Solution for the Case p = 1 

We here consider exactly the same problem as in Subsection 16.1.21 i.e., 
our parabolic model problem is posed in the 3d space-time domain Q = 
(0,1)^. Eor simplicity, we here consider only the case p = 1 that turns 
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Fig. 7. Solution contours for Example I6.2..SI (left) and the control points 
(right) are given by Pii,* 2 ,i 3 = { (0,0,0), (1,0,0), (0,1,0), (1,1,0), 
(0.25,0,0.50), (0.75,0,0.5), (0.25,1,0.5), (0.75,1,0.5), (0,0,1), (1,0,1), 
(0,1,1), (1,1,1)} for ii = 1, 2, i 2 = 1,2 and h = 1, 2, 3. 


p = 3 

p = 4 

Dofs 

II/i,m 

Rate 

Dofs 


Rate 

64 

2.41099e-01 

0 

125 

9.93223e-03 

0 

125 

2.65597e-02 

3.18 

216 

7.60959e-03 

0.38 

343 

4.25800e-03 

2.64 

512 

4.70558e-04 

4.02 

1331 

4.58251e-04 

3.21 

1728 

2.72696e-05 

4.11 

6859 

5.46583e-05 

3.06 

8000 

1.60187e-06 

4.09 

42875 

6.73180e-06 

3.02 

42875 

9.81655e-08 

4.03 


Table 18. Errors and rates wrt 
p = 3 and p = 4. 


h.m 


for Example 16.2 


3] using degrees 


P = 1 

p = 2 

Dofs 

I|m - Uh\\L2(Q) 

Rate 

Dofs 

Ih - 'Uh\\L2(Q) 

Rate 

8 

3.37961e-01 

0 

27 

2.67611e-02 

0 

27 

7.67535e-02 

2.14 

64 

2.66302e-02 

0.01 

125 

1.95993e-02 

1.96 

216 

2.33938e-03 

3.51 

729 

4.86526e-03 

2.01 

1000 

2.63360e-04 

3.15 

4913 

1.25329e-03 

1.95 

5832 

3.20369e-05 

3.04 

35937 

3.61179e-04 

1.79 

39304 

3.97751e-06 

3.00 


Table 19. L 2 errors and rates for Examp 


degrees p 


p = 2. 


1 and 
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p = 3 

p = 4 

Dots 

M - Uh\\L2(Q) 

Rate 

Dofs 

\\u-Uh\\L2(Q) 

Rate 

64 

2.65683e-01 

0 

125 

6.83584e-04 

0 

125 

2.26323e-03 

3.55 

216 

6.75251e-04 

0.02 

343 

3.36653e-04 

2.74 

512 

3.96777e-05 

4.09 

1331 

1.70973e-05 

4.30 

1728 

1.60770e-06 

4.63 

6859 

9.93566e-07 

4.10 

8000 

5.11255e-08 

4.97 

42875 

6.12843e-08 

4.02 

46656 

1.61035e-09 

4.99 


Table 20. L 2 errors and rates for Examp 
p = 4. 


e 16.2.31 using degrees p = 3 and 


out to be the finite element case with trilinear hexahedral finite elements. 
First, we decompose the space-time mesh into several subdomains and 
assemble the arising linear systems of Subsection 16.1.21 in parallel. For 
example, in Figure [SI the space-time decomposition with 64 subdomains 
is shown. Afterwards, we solve these linear systems also in parallel with 
the GMRES method, where we apply the AMG library hypre as a pre¬ 
conditioner. For the stopping criteria, we use the relative residual error 
reduction by 10“^^. In Table EH we show the iteration numbers and the 
solving times for different uniform refinement levels, where we increase 
the number of cores for larger problems. We observe that the iteration 
numbers are slightly increasing, but the performance of this solver is still 
very good, even if this solver was not constructed for space-time problems. 
This example was computed on the supercomputer Vulcan BlueGene/Q 
in Livermore, Galifornia U.S.A by using the finite element library MFEM. 


Table 


Dofs 

||w - MJilUaCQ) 

Rate 

iter 

time [s] 

cores 

8 

3.65528e-01 

- 

1 

0.01 

1 

27 

9.39008e-02 

1.961 

2 

0.01 

1 

125 

2.32674e-02 

2.013 

6 

0.01 

1 

729 

5.75635e-03 

2.015 

15 

0.07 

64 

4 913 

1.43198e-03 

2.007 

16 

0.14 

64 

35 937 

3.57217e-04 

2.003 

19 

0.40 

64 

274 625 

8.92171e-05 

2.001 

24 

1.04 

1 024 

2 146 689 

2.22941e-05 

2.001 

29 

3.65 

1 024 

16 974 593 

5.57231e-06 

2.000 

36 

21.40 

1 024 

135 005 697 

1.39293e-06 

2.000 

50 

36.26 

8 192 

1 076 890 625 

3.48206e-07 

2.000 

63 

156.50 

16 384 


21. Solver performance for Example 16.1.21 and p = 


1 . 
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Fig. 8. Space-time decomposition with 64 subdomains. 


7 Conclusions 

We have introduced the Space-Time IgA for parabolic evolution problems. 
We have presented a-priori error estimates and numerical examples in the 
space-time domain Q = {{x,t) £ M'^+^ : x G G (0,T)} for both 

hxed spatial domains 17 C and moving spatial domains 17 (t) C 
t G [0,T]. Our numerical experiments have been preformed on a sequence 
of refined meshes and for the polynomial degrees p = 1,2,3,4 in the 
cases d = 1 as well as d = 2. In the case of smooth solutions, we have 
nicely observed the full asymptotical convergence rates. For simplicity, 
we restricted our-self to the single-patch case. However, it is possible 
to generalize the results to the conform multi-patch case. Moreover, the 
combination of the results of this paper with the results of the papers [33] 
and |39] allows us to analyze space-time multi-patch dG IgA schemes in a 
similar way. The overall efficiency of the space-time IgA heavily depends 
on the availability of fast parallel solvers. At the first glance, the solution 
of one large space-time system of linear algebraic equations instead of 
many smaller systems in traditional time-stepping methods seems to be 
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a big disadvantage of space-time IgA, but on parallel computers with 
many cores this is a big advantage that allows us to overcome the curse 
of sequentiality as Example 16.31 shows. Another advantage consists in 
the elegant treatment of moving domains or interfaces. And last, but 
not least the space-time adaptivity or, more precisely, the possibility to 
perform a free adaptivity in Q without separating the time from the space 
{t is just another variable Xd+i) opens new horizons in developing highly 
efficient parallel adaptive space-time IgA methods for parabolic as well 
as hyperbolic problems, using T-splines [l5] or THB-splines |16] for local 
refinement. 
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